{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "progressive-travel",
   "metadata": {},
   "source": [
    "# Limits to Growth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "black-toolbox",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "earlier-pride",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bound-nature",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "found-pledge",
   "metadata": {},
   "source": [
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "Click here to access the notebooks: <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "general-noise",
   "metadata": {
    "tags": []
   },
   "source": [
    "Here's the data from the previous chapter again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "affiliated-eleven",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/data/World_population_estimates.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "recent-trouble",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pandas import read_html\n",
    "\n",
    "filename = 'World_population_estimates.html'\n",
    "tables = read_html(filename, header=0, index_col=0, decimal='M')\n",
    "table2 = tables[2]\n",
    "table2.columns = ['census', 'prb', 'un', 'maddison', \n",
    "                  'hyde', 'tanton', 'biraben', 'mj', \n",
    "                  'thomlinson', 'durand', 'clark']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "western-blowing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "un = table2.un / 1e9\n",
    "census = table2.census / 1e9"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "occasional-kitchen",
   "metadata": {
    "tags": []
   },
   "source": [
    "And here are the functions from the previous chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "simple-coupon",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'chap06.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "monetary-profile",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from chap06 import run_simulation\n",
    "\n",
    "def plot_estimates():\n",
    "    census.plot(style=':', label='US Census')\n",
    "    un.plot(style='--', label='UN DESA')\n",
    "    decorate(xlabel='Year', \n",
    "             ylabel='World population (billions)') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "damaged-reservation",
   "metadata": {},
   "source": [
    "In the previous chapter we developed a population model where net growth during each time step is proportional to the current population. This model seems more realistic than the constant growth model, but it does not fit the data as well.\n",
    "\n",
    "There are a few things we could try to improve the model:\n",
    "\n",
    "-   Maybe net growth depends on the current population, but the\n",
    "    relationship is quadratic, not linear.\n",
    "\n",
    "-   Maybe the net growth rate varies over time.\n",
    "\n",
    "In this chapter, we'll explore the first option.\n",
    "In the exercises, you will have a chance to try the second. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "assigned-slovakia",
   "metadata": {},
   "source": [
    "## Quadratic Growth\n",
    "\n",
    "It makes sense that net growth should depend on the current population, but maybe it's not a linear relationship, like this:\n",
    "\n",
    "```\n",
    "net_growth = system.alpha * pop\n",
    "```\n",
    "\n",
    "Maybe it's a quadratic relationship, like this:\n",
    "\n",
    "```\n",
    "net_growth = system.alpha * pop + system.beta * pop**2\n",
    "```\n",
    "\n",
    "We can test that conjecture with a new update function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "beginning-belly",
   "metadata": {},
   "outputs": [],
   "source": [
    "def growth_func_quad(t, pop, system):\n",
    "    return system.alpha * pop + system.beta * pop**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "initial-factory",
   "metadata": {},
   "source": [
    "Here's the `System` object we'll use, initialized with `t_0`, `p_0`, and `t_end`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "listed-florence",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = census.index[0]\n",
    "p_0 = census[t_0]\n",
    "t_end = census.index[-1]\n",
    "\n",
    "system = System(t_0=t_0,\n",
    "                p_0=p_0,\n",
    "                t_end=t_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amber-context",
   "metadata": {},
   "source": [
    "Now we have to add the parameters `alpha` and `beta` .\n",
    "I chose the following values by trial and error; we'll see better ways to do it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "signed-impossible",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.alpha = 25 / 1000\n",
    "system.beta = -1.8 / 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "confidential-retreat",
   "metadata": {},
   "source": [
    "And here's how we run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "italian-converter",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, growth_func_quad)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "forbidden-brisbane",
   "metadata": {},
   "source": [
    "Here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "simplified-sight",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.plot(color='gray', label='model')\n",
    "plot_estimates()\n",
    "decorate(title='Quadratic growth model')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "primary-ending",
   "metadata": {},
   "source": [
    "The model fits the data well over the whole range, with just a bit of space between them in the 1960s.\n",
    "\n",
    "It is not entirely surprising that the quadratic model fits better than the\n",
    "constant and proportional models, because it has two parameters we can\n",
    "choose, where the other models have only one. In general, the more\n",
    "parameters you have to play with, the better you should expect the model\n",
    "to fit.\n",
    "\n",
    "But fitting the data is not the only reason to think the quadratic model\n",
    "might be a good choice. It also makes sense; that is, there is a\n",
    "legitimate reason to expect the relationship between growth and\n",
    "population to have this form.\n",
    "\n",
    "To understand it, let's look at net growth as a function of population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sunset-underground",
   "metadata": {},
   "source": [
    "## Net Growth\n",
    "\n",
    "Let's plot the relationship between growth and population in the quadratic model.\n",
    "I'll use `linspace` to make an array of 101 populations from 0 to 15 billion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "neural-guinea",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace\n",
    "\n",
    "pop_array = linspace(0, 15, 101)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "heated-selling",
   "metadata": {},
   "source": [
    "Now I'll use the quadratic model to compute net growth for each population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "animal-spoke",
   "metadata": {},
   "outputs": [],
   "source": [
    "growth_array = (system.alpha * pop_array + \n",
    "                system.beta * pop_array**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "engaging-parade",
   "metadata": {},
   "source": [
    "To plot growth rate versus population, we'll use the `plot` function from Matplotlib.\n",
    "First we have to import it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "informed-three",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "retained-deployment",
   "metadata": {},
   "source": [
    "Now we can use it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "unexpected-nigeria",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(pop_array, growth_array, label='net growth', color='C2')\n",
    "\n",
    "decorate(xlabel='Population (billions)',\n",
    "         ylabel='Net growth (billions)',\n",
    "         title='Net growth vs. population')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "precise-finish",
   "metadata": {},
   "source": [
    "Note that the x-axis is not time, as in the previous figures, but population. We can divide this curve into four kinds of behavior:\n",
    "\n",
    "-   When the population is less than 3 billion, net growth is\n",
    "    proportional to population, as in the proportional model. In this\n",
    "    range, the population grows slowly because the population is small.\n",
    "\n",
    "-   Between 3 billion and 10 billion, the population grows quickly\n",
    "    because there are a lot of people.\n",
    "\n",
    "-   Above 10 billion, population grows more slowly; this behavior models\n",
    "    the effect of resource limitations that decrease birth rates or\n",
    "    increase death rates.\n",
    "\n",
    "-   Above 14 billion, resources are so limited that the death rate\n",
    "    exceeds the birth rate and net growth becomes negative.\n",
    "\n",
    "Just below 14 billion, there is a point where net growth is 0, which\n",
    "means that the population does not change. At this point, the birth and death rates are equal, so the population is in *equilibrium*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "angry-voice",
   "metadata": {},
   "source": [
    "## Finding Equilibrium\n",
    "\n",
    "The equilibrium point is the population, $p$, where net population growth, $\\Delta p$, is 0.\n",
    "We can compute it by finding the roots, or zeros, of this equation: \n",
    "\n",
    "$$\\Delta p = \\alpha p + \\beta p^2$$ \n",
    "\n",
    "where $\\alpha$ and $\\beta$ are the parameters of the model. \n",
    "If we rewrite the right-hand side like this: \n",
    "\n",
    "$$\\Delta p = p (\\alpha + \\beta p)$$ \n",
    "\n",
    "we can see that net growth is $0$ when $p=0$ or $p=-\\alpha/\\beta$.\n",
    "So we can compute the (non-zero) equilibrium point like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ordinary-honolulu",
   "metadata": {},
   "outputs": [],
   "source": [
    "-system.alpha / system.beta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adaptive-pharmacy",
   "metadata": {},
   "source": [
    "With these parameters, net growth is 0 when the population is about 13.9 billion\n",
    "(the result is positive because `beta` is negative).\n",
    "\n",
    "In the context of population modeling, the quadratic model is more\n",
    "conventionally written like this: \n",
    "\n",
    "$$\\Delta p = r p (1 - p / K)$$ \n",
    "\n",
    "This is the same model; it's just a different way to *parameterize* it. Given $\\alpha$ and $\\beta$, we can compute $r=\\alpha$ and $K=-\\alpha/\\beta$.\n",
    "\n",
    "In this version, it is easier to interpret the parameters: $r$ is the\n",
    "unconstrained growth rate, observed when $p$ is small, and $K$ is the\n",
    "equilibrium point. \n",
    "$K$ is also called the *carrying capacity*, since it indicates the maximum population the environment can sustain."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continental-image",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this chapter we implemented a quadratic growth model where net growth depends on the current population and the population squared.\n",
    "This model fits the data well, and we saw one reason why: it is based on the assumption that there is a limit to the number of people the Earth can support.\n",
    "\n",
    "In the next chapter we'll use the models we have developed to generate\n",
    "predictions.\n",
    "But first, I want to warn you about a few things that can go wrong when you write functions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eligible-pride",
   "metadata": {},
   "source": [
    "## Dysfunctions\n",
    "\n",
    "When people learn about functions, there are a few things they often\n",
    "find confusing. In this section I'll present and explain some common\n",
    "problems.\n",
    "\n",
    "As an example, suppose you want a function that takes a\n",
    "`System` object, with variables `alpha` and `beta`, and computes the\n",
    "carrying capacity, `-alpha/beta`. \n",
    "Here's a good solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "realistic-opinion",
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olive-information",
   "metadata": {},
   "source": [
    "Now let's see all the ways that can go wrong."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prostate-motorcycle",
   "metadata": {},
   "source": [
    "*Dysfunction #1:* Not using parameters. In the following version, the function doesn't take any parameters; when `sys1` appears inside the function, it refers to the object we create outside the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "marine-entry",
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity():\n",
    "    K = -sys1.alpha / sys1.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity()\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dated-invalid",
   "metadata": {},
   "source": [
    "This version works, but it is not as versatile as it could be.\n",
    "If there are several `System` objects, this function can work with only one of them, and only if it is named `sys1`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "meaningful-louisiana",
   "metadata": {},
   "source": [
    "*Dysfunction #2:* Clobbering the parameters. When people first learn\n",
    "about parameters, they often write functions like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "moving-brazil",
   "metadata": {},
   "outputs": [],
   "source": [
    "# WRONG\n",
    "def carrying_capacity(system):\n",
    "    system = System(alpha=0.025, beta=-0.0018)\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.03, beta=-0.002)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dietary-spectacular",
   "metadata": {},
   "source": [
    "In this example, we have a `System` object named `sys1` that gets passed\n",
    "as an argument to `carrying_capacity`. But when the function runs, it\n",
    "ignores the argument and immediately replaces it with a new `System`\n",
    "object. As a result, this function always returns the same value, no\n",
    "matter what argument is passed.\n",
    "\n",
    "When you write a function, you generally don't know what the values of\n",
    "the parameters will be. Your job is to write a function that works for\n",
    "any valid values. If you assign your own values to the parameters, you\n",
    "defeat the whole purpose of functions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "present-estonia",
   "metadata": {},
   "source": [
    "*Dysfunction #3:* No return value. Here's a version that computes the value of `K` but doesn't return it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "sacred-physiology",
   "metadata": {},
   "outputs": [],
   "source": [
    "# WRONG\n",
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "technological-incentive",
   "metadata": {},
   "source": [
    "A function that doesn't have a return statement actually returns a special value called `None`, so in this example the value of `pop` is `None`. If you are debugging a program and find that the value of a variable is `None` when it shouldn't be, a function without a return statement is a likely cause."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "received-firewall",
   "metadata": {},
   "source": [
    "*Dysfunction #4:* Ignoring the return value. Finally, here's a version where the function is correct, but the way it's used is not.\n",
    "\n",
    "```\n",
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "carrying_capacity(sys1)   # WRONG\n",
    "print(K)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "approximate-straight",
   "metadata": {},
   "source": [
    "In this example, `carrying_capacity` runs and returns `K`, but the\n",
    "return value doesn't get displayed or assigned to a variable.\n",
    "If we try to print `K`, we get a `NameError`, because `K` only exists inside the function.\n",
    "\n",
    "When you call a function that returns a value, you should do something\n",
    "with the result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "liable-mixture",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worst-builder",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    " In a previous section, we saw a different way to parameterize the quadratic model:\n",
    "\n",
    "$$ \\Delta p = r p (1 - p / K) $$\n",
    "\n",
    "where $r=\\alpha$ and $K=-\\alpha/\\beta$.  \n",
    "\n",
    "Write a version of `growth_func` that implements this version of the model.  Test it by computing the values of `r` and `K` that correspond to `alpha=0.025` and `beta=-0.0018`, and confirm that you get the same results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "stretch-check",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "tender-treat",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "passive-certificate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "understood-cancer",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "  What happens if we start with an initial population above the carrying capacity, like 20 billion?  Run the model with initial populations between 1 and 20 billion, and plot the results on the same axes.\n",
    "\n",
    "Hint: If there are too many labels in the legend, you can plot results like this:\n",
    "\n",
    "```\n",
    "    results.plot(label='_nolegend')\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "agricultural-burke",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "colored-globe",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
